Hydrodynamic fluctuations in confined emulsions 
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When an ensemble of particles interact hydrodynamically, they generically display large-scale 
transient structures such as swirls in sedimenting particles [l] , or colloidal strings in sheared suspen- 
sions [2J. Understanding these nonequilibrium fluctuations is a very difficult problem, yet they are 
of great importance for a wide range of processes including pigment deposition in cosmetic or paint 
films, the transport of microfluidic droplets, ... All these samples concern rigidly confined fluids, 
which we consider in this paper. We address the collective dynamics of non-Brownian droplets cruis- 
ing in a shallow microchannel. We provide a comprehensive characterization of their spatiotemporal 
density fluctuations. We show that density excitations freely propagate at all scales, and in all direc- 
tions even though the particles are neither affected by potential forces nor by inertia. We introduce 
a theory which quantitatively accounts for our experimental findings. By doing so we demonstrate 
that the fluctuation spectrum of this nonequilibrium system is shaped by the combination of truly 
long-range hydrodynamic interactions and local collisions. 



Understanding the collective dynamics of non- 
brownian particles in viscous fluids is a long-standing 
challenge in fluid mechanics. A typical example is the 
still unsolved problem of sedimentation in a quiescent 
fluid. Rather than falling along straight lines, as an 
isolated particle would do, sedimenting particles expe- 
rience additional swirling motion correlated over large 
yet finite distances. The physical origin of those dynam- 
ical structures has been under debate for more than 30 
years [31 0]. The conceptual complexity of this collec- 
tive dynamics contrasts with the formal simplicity of the 
(linear) Stokes equation that rules low-Reynolds-number 
flows. Immersed bodies generically affect both the mo- 
mentum and the mass transfers of the fluid, even when 
not driven by external fields. As a result, long-range 
effective interactions arise between the particles due to 
the interplay between the local velocity of the fluid and 
the motion of the particles. The particle transport prob- 
lem thus corresponds to a many-body dynamical system, 
that is intrinsically non-linear. The hydrodynamic inter- 
actions actually vanish only for uniform flow fields, for 
which the particles would be all advected at the same 
speed as the fluid, irrespective of their spatial distribu- 
tion. Such a condition is never achieved when the fluid is 
rigidly confined. For instance, when particles are driven 
in thin channels or in (semi)rigid films, the momentum 
exchange with the bounding walls causes strong distor- 
tions of the flow field, thereby inducing effective inter- 
actions between the particles [SHE]- As it turns out, the 
transport of particle-ladden fluid in thin films is involved 
in a number of industrial and natural processes, including 
protein motion in lipid membranes [5], bacteria swarm- 
ing [TO] j colloid deposition on solid surfaces [TTJ [T2] , 
and droplet-based microfluidics [T3]. Understanding the 



particle transport in rigidly confined films is a first step 
toward the description of particle traffic in more complex 
geometries such as ordered, or random porous networks. 
The pioneering experiments by Rouyer et al. revealed 



A.iU 


■ i 

'________« 






X 



FIG. 1. A- Picture of the microfluidic setup. It is composed 
of: (i) a conventional flow-focusing junction, (ii) two dilution 
channels, and (iii) two identical wide observation channels. 
During the experiment one of these two 5-cm long channels is 
continuously fed with monodisperse droplets having a diam- 
eter comparable to the channel height. Scale bar: 5 cm. B- 
Close up on the main channel. Droplet diameter: 33.4 /im. 
Scale bar: 5 mm. C- Typical snapshot of a movie used to 
track the droplet positions. The field of view corresponds 
to the dashed rectangle in B. The black arrows indicate the 
direction of the flow. Scale bar: 500 /im. 
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that an ensemble of confined particles advected by a uni- 
form viscous flow displays short-range structural corre- 
lations akin to the one found in simple liquids at the 
molecular scale [TJ]. More recently, Beatus and cowork- 
ers exploited a simple microfluidic experiment to probe 
the propagation of density heterogeneities in bidimen- 
sional emulsions [T5]. They showed that the dynamic 
creation and destruction of droplet clusters also display 
a temporal coherence. Focusing on a semi-local quantity, 
the area fraction averaged over the channel width, they 
demonstrated that the droplet velocity is an increasing 
function of the local density, which is the elementary in- 
gredient to give rise to Burgers (nonlinear) shock waves. 
The same phenomenology was also found in focused par- 
ticle streams [15]. However, this elegant observation does 
not account for the complexity of the spatiotemporal fluc- 
tuations generically observed at all scales in rigidly con- 
fined particle-laden fluids. 

Here we combine experimental and theoretical results 
to shed light on the collective dynamics of particles cruis- 
ing through large-scale microfluidic channels. We first 
provide a comprehensive characterization of the density 
fluctuations observed in our experiments, and unveil their 
propagative nature at all scales, and in all directions. 
Then, we quantitatively demonstrate how the fluctua- 
tion spectrum of this many-body non-equilibrium system 
is shaped by the combination of truly long-range hydro- 
dynamic interactions and local collisions. 

EXPERIMENTAL RESULTS 

Briefly, our experimental system consists of a monodis- 
perse emulsion flowing in a shallow microchannel. The 
length and width of the channel, L x W — 5 cm x 5 mm, 
are much larger than its height, h = 27 ± 0.1 pm, 
which compares with the droplet diameter, see Fig. 1. 
The emulsion is therefore confined in a quasi-2D geom- 
etry. The droplets are formed at a conventional flow- 
focusing junction followed by a dilution module. The 
fluid flow-rates are imposed by high-precision syringe 
pumps. Etched-glass microchips ensure that the chan- 
nel dimensions are unaffected by the flow conditions. In 
addition, the geometry of the junction, and the range 
of flow rates, are chosen so-that the formation of the 
droplet was unaffected by the dilution flow. By doing so, 
we accurately control both the droplet radius, Rd, and 
the average area fraction, </>, occupied by the emulsion. 
We report here results obtained for = 16.7 ± 0.3 pm 
(R d /h = 0.62), and 0.21 < (j> < 0.56. Varying the droplet 
sizes up to i?d ~ 2h does not qualitatively change our 
measurements. The droplets are visualized using fluores- 
cence imaging. For each experiment we tracked ~ 10 5 
particle trajectories in a region close to the center of the 
main channel, see Fig. IB and 1C. More details about 
the microfluidic and the imaging setup are provided in 



the Materials and Methods section. 

In the absence of droplets, the fluid flow would be uni- 
form along the x-direction in the observation region. This 
is evidenced by the linear trajectories followed by iso- 
lated droplets cruising along the channel. Conversely, 
even at the smallest surface fraction, when an emulsion 
flows, the droplets undergo large transverse and longitu- 
dinal fluctuations in their motion, as shown in the sup- 
plementary movie (Video SI). These fluctuations lead to 
the formation of particle clusters at all scales. Interest- 
ingly, these clusters are clearly seen to travel at a speed 
that is different from the mean droplet velocity: trans- 
verse clusters are faster, whereas the longitudinal ones 
are slower. However, these clusters are transient struc- 
tures, they form and break apart in a continuous fashion. 
Our purpose is to elucidate the physical mechanisms re- 
sponsible for this complex and fluctuating dynamics. To 
quantify the spatiotemporal fluctuations of the droplet 
density field p{r,t), where r = (x,y), we measure its 
power spectrum, that is the dynamic structure factor. 
Introducing the Fourier transform of the local density: 
Pq,w' = 2^f / /°( r > *) e ^ q r ~" '^dr dt, the power spectrum 
is defined as |/5q, W '| 2 , where p(r,t) = p(r,t) — (p(r,t)). 
Practically, p is computed from the particle positions as 
p(r, t) = J^i G( r ~ r i(t)), where rj(i) is the position of the 
i th droplet, and where Q is a Gaussian shape function of 
width i?d/15, see Materials and Methods. 

In Fig 2A, we show a slice of a typical power spec- 
trum in the (u/ ,q x ) plane. This example corresponds 
to <f> = 0.39, and to q y Rd = 0.2. Several important 
comments are in order: (i) The power spectrum is lo- 
calized in the Fourier space, which is the hallmark of 
propagative dynamics for the density fluctuations, as first 
noted in [TS] for the specific case of the y-averaged den- 
sity mode {q y = 0). We stress that compression modes 
propagate even though the droplets do not interact via 
potential forces, and even though their inertia is negligi- 
ble compared to the viscous friction at this scale. These 
"sound" modes originate only from the hydrodynamic 
coupling between the advected particles, (ii) The curve 
on which the spectrum is peaked corresponds to the dis- 
persion curve of the density waves. We emphasize that it 
deviates markedly from a straight line at moderate wave- 
lengths. The hydrodynamic interactions do not merely 
renormalizc the mean advection speed and cause the den- 
sity fluctuations to propagate in a dispersive fashion, (iii) 
The global shape of the spectrum is conserved for every 
area fraction, and more surprisingly for every wave vec- 
tor q y provided that the wavelength remains larger than 
the particle size (see below). 

In all that follows, we discard the trivial non-dispersive 
contribution due to the advection at the mean droplet ve- 
locity (vd). To do so, we focus on the density fluctuations 
in the frame moving at (v<j), and introduce the reduced 
pulsation u = u/ — (vd)q x - Experiments done at differ- 
ent area fractions, and thus at different continuous phase 
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FIG. 2. A- Grayscale power spectrum of the density fluctuations plotted in the (q x ,u)') plane for q y — 0.2/Rd- Area fraction 
4> — 0.39, mean fluid velocity: vf — lmm/s. Solid line: theoretical prediction for the location of the dispersion curve. 
B- Experimental dispersion curve in the frame moving at (vd), for 4> = 0.39. uj is plotted versus q x , and q y . We recall 
that units are chosen so that Rd = 1, and Wp = 1. The dotted line indicates the q y value corresponding to the power 
spectrum shown in A. C- Renormalized dispersion relations in the moving frame, for q y — 0. Circles: experimental data (<f> = 
0.21, 0.25, 0.30, 0.31, 0.33, 0.35, 0.38, 0.39, 0.41, 0.45, 0.49, 0.56). Solid line: theoretical prediction, Eq. |8j with no adjustable 
parameter. D- Variations of the (rescaled) mean droplet velocity as a function of the area fraction. Circles: experimental data. 
Solid line: best linear fit. The error bars account for statistical fluctuations, and correspond to the standard deviation. E- w gx 
plotted versus q v at q x = 0. Circles: experimental data for <j) = 0.56. Solid line: Theoretical prediction with no adjustable 
parameter deduced from Eq. [8] The error bars correspond to a 95% confidence interval in the measurement of v gx from the 
slope of the dispersion curve. F- Theoretical prediction for the variations of the group velocity, u gx with the wave vector 
components. Dotted line: q x — 0. The variations of « gx along this direction are shown in E. 



velocities, are compared by normalizing the pulsations, 
and the wave vectors by v-p/Rd, and R^ 1 respectively. 
Fig 2B shows a typical dispersion relation: oj — u(q x , q y ), 
obtained for <fi = 0.39. The spectrum is symmetric along 
the q y direction as expected from the symmetry of the 
system. More importantly, density fluctuations propa- 
gate in all directions except in the one strictly transverse 
to the flow (q x =0). In addition, the dispersion curve 
displays an axial symmetry with respect to the q y -axis. 
It is worth noting that the sign of the associated phase 
velocity changes as q x increases. The long wavelength 
excitations propagate downstream, while the short wave- 
length excitations propagate upstream. 

In Fig. 2C, we show that once renormalized by </>, 
the dispersion relations corresponding to 12 different area 
fractions collapse on a single master curve. This notice- 



able collapse is not specific to the purely longitudinal 
waves and occurs for all the possible q y values. This 
systematic rescaling demonstrate that a unique set of 
physical mechanisms dictates the collective motion of the 
droplets, at all scales, regardless of the droplet density. 
We now establish a theoretical model which quantita- 
tively accounts for these robust experimental findings. 



THEORY AND DISCUSSION 

We did not observe any deformation of the droplets in 
our experiments. Therefore, the instantaneous configura- 
tion of the emulsion is fully determined by the positions 
of N identical axisymmetric particles: rj(i), i — 1 . . . N. 
The dynamics of an isolated particle has proven to be 
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correctly captured by a constant mobility coefficient, p, 
defined as r,-(t) = //v(rj,i) where v(r, t) is the in-plane 
fluid velocity field averaged over the channel height in 
the absence of the particle i [T5] . In our quasi-2D geome- 
try, the fluid flow is a potential flow and derives from the 
local pressure field: v = — GVP where G = h 2 /12i], rj 
being the viscosity of the aqueous phase. v(r, t) is then 
fully determined when considering the incompressibility 
condition, and the no-flux boundary conditions through 
the sidewalls of the channel as well. In a particle-free 
channel, the velocity field would be uniform, v = upx. 
The particles are not passive tracers (p < 1), there- 
fore their relative motion with respect to the fluid re- 
sults in a dipolar disturbance of the surrounding flow, 
as shown in [5J [B] . The potential dipolar perturbation, 
v dlp (r, Ti(i)), induced at the position r by a particle lo- 
cated at Ti (t) is defined by the modified incompressibility 
relation: 



dip 



(v,r i (t))=ad x S(r-v l (t)) 



(1) 



where a is the dipole strength (er > 0). To establish the 
equations of motion of the N particles, we now assume 
the dipolar disturbances to be pairwise additive. This 
yields: 



r<(t) = pvpi + pj^^' 



M<),r,(i)) 



(2) 



We now show how to move from these N coupled equa- 
tions to an hydrodynamic description for the particle 
density field p(r,t). We follow a well established kinetic 
theory framework, see e.g. [Ill [18]. We first introduce the 
iV-point distribution function, p( N \ri, Tff, t), to find 
particles at r lv .., rjy at the time t. In a stationary state 



di 



0. We use this relation and Eq. [ljto derive the 



Fokkcr-Planck equation obeyed by p^ : 



JY 



(N) 



(3) 



5> di P(r^),r^))p 



We now exploit two assumptions: the hydrodynamic 
interactions are assumed to be pairwise additive, and the 
particles are all identical. Therefore we can establish a 
mass conservation equation that couples the local density 
p(r,t) to the two-point distribution function p^> (r, r', t) 
only. It is obtained by integrating Eq. [3] over N — 1 
particle positions; 

d t p(r, t) + pv F d x p(r, t) (4) 

+ pV ■ J dr'v dip (r,r')p (2) (r,r',t) =0 

where, p(r,t) = {N l iy J p (JV) (r,r 2 , ...,r N ,t)dr 2 ■ •■dr JV , 
p( 2 )(r,r',i) = w h W fpW(r,r>,r 3 ,...,r N ,t)dr 3 ...dr N . 



We now assume that the particle positions decorrelate 
over a distance as small as one particle diameter. In ad- 
dition to this mean-field approximation, we also explic- 
itly account for the steric repulsion between the particles. 
To do so, we postulate the following closure relation for 
Eq.gj 



P (2) (r, 





p(r)p(r') 
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<2i? d 
> 2R d 



(5) 



where i?d is the radius of a particle. Eqs. [4] and [5j de- 
fine the equations of motion for the particle-density field. 
In principle, the effective extent of the excluded volume 
could be larger that the particle radius due to short-range 
intermolecular repulsions, and lubrication forces. How- 
ever no measurable difference with the actual droplet ra- 
dius could be observed in our experiments. 

We now focus on the dynamics of small density fluctu- 
ations, p(r,t), around an homogeneous state: p(r,t) = 
p(r,t) - po, where p = (p(r,t)) = 4>/{ttRI). As 
done in our experiments, we work in the frame mov- 
ing at the mean droplet velocity (v d ) = pv F x + 
Wo J\r-r'\>2R d v dlp (r , r') dr' . At leading order in p, 
Eqs. [4]ancr[5] then reduce to: 



ftp(r,t)+V-j(r,t) = 



(6) 



where both the hydrodynamic interactions (long- 
range) and the contact interactions (short-range) 
are captured by the current functional j(r, t) = 
Wo I\r-r'\>2R d v dip (r, r')p(r', f) dr'. Using Eq. 1 and fo- 
cusing on particles far from the sidewalls, we can show 
that V • j is actually a local quantity: 



V-j(r,i) 
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p(r-2R d r')cos6'd6' (7) 



Since i? d "C W we have used the expression of the 
dipolar perturbation corresponding to an unbounded do- 
main: v dip (r,r + 2i? d r') • f = -(a cos ff)/%nR\, where 
f ' = cos 9'x + sin 9'y [H] . We now look for plane wave 
solutions p(r, t) — p q exp(iujt — iq • r) of Eq. |6] After 
some elementary algebra we infer their dispersion rela- 
tion, which is our main theoretical result: 



oj = {pap )q x 



Ji(2gR d ) 
2qR d 



(8) 



where J± is the first Bessel function. As uj is a real quan- 
tity, the above equation implies that density waves freely 
propagate in the channel in qualitative agreement with 
our experimental observations. It is worth noting that 
since V • j is a local quantity, the form of the dispersion re- 
lation is generic, and does not depend on the channel size, 
and geometry. In addition, we point that uj scales as po, 
which explains the collapse of the normalized dispersion 
relations on a single master curve over the entire range of 
wave vectors, Figure [2p. We now move to a quantitative 
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comparison between our theoretical prediction and our 
experimental measurements. Eq. [8] is fully determined 
by two physical parameters: the droplet radius Rd, and 
pcrpo that quantifies the strength of the hydrodynamic 
couplings. To determine this latter parameter, we exploit 
another specific feature of the long-range hydrodynamic 
interactions. In an isotropic and homogeneous system, 
due to their dipolar symmetry, the sum of all the hy- 
drodynamic interactions would leave the mean droplet 
velocity unchanged. However, in anisotropic channel ge- 
ometries, (vd) increases linearly with the mean density 
irrespective of the channel size [19 . At th order in p, 
(v c i) = pvpi + 5(/iO"po)x. Importantly this relation pro- 
vides a direct means to measure independently the last 
unknown parameter of our theory. Figure [2p clearly 
shows that the measured mean droplet speed increases 
with po in an affinc manner. The strength of the hydrody- 
namic coupling (papo) is readily infered from a linear fit, 
Figure [2jD. We superimposed our theoretical predictions 
for the dispersion relation, Eq. [8] both in the laboratory 
frame and in the frame moving at (vd) in Figures [2|\ 
and [2p. We find that the agreement between the the- 
ory and the experiments is excellent over a wide range 
of wave vectors, and of area fractions. Without any free 
fitting parameters, our model quantitatively captures the 
dispersive nature of the density fluctuations observed in 
the flowing emulsions. 

To gain more physical insight into the propagation of 
the density waves, it is worth looking at the small-g ex- 
pansion of Eq.M u> = \po-p q x [l— ^(qR d ) 2 ]+0 ((qRd) 4 )- 
At leading order, this relation is non-dispersive (linear) 
whatever the direction of propagation. The phase ve- 
locity scales linearly with the magnitude of the dipolar 
coupling a. In addition, it does not depend explicitly 
on Rd, which implies that the small-g excitations prop- 
agate only due to the long-range hydrodynamic interac- 
tions between the particles. Conversely, the dispersive 
term in w(q) explicitly depends on the particle radius. 
At high q, the propagation of the density waves is set by 
the combination of the excluded volume interactions and 
the angular symmetry of the hydrodynamic couplings. 

We close this section by recalling that one of the most 
striking feature observed in the flowing emulsions is the 
propagation of vertical density bands which propagate 
significantly faster than the mean droplet flow, see sup- 
plemental material (Video SI). An homogeneous verti- 
cal band spanning the entire width of the channel corre- 
sponds to the linear superposition of plane waves associ- 
ated with q y — 0, and with q x s distributed around q x = 0. 
In the frame moving at (vd), their speed is given by the 
x-component of the group velocity v g3 .(q x ,q y ) = dw/dq x 
evaluated at q = 0. In Figure (2}5, we plot the experi- 
mental values of v ga .(0, q y ), which we measured from the 
slope at the origin of the dispersion curves (as the ones 
shown in Figure [2p). Again the agreement with the the- 
oretical curve deduced from Eq. [8] is excellent. This plot 



reveals that the density bands extended across the entire 
channel width are the fastest and propagate at velocities 
1.5 higher that the mean droplet flow, thereby making 
them highly visible on the experimental movies. 

To go beyond this observation, we use Eq. U to plot 
the magnitude of v ga .(q x ,qy) for all qs, see Fig. [2j?. v gx 
displays non-monotonic variations with both q x , and q y 
and changes its sign at high qs. In the frame moving at 
(vd), and for qRd -C 1, v Sx is positive. The wave packets 
propagate faster than the mean flow due to the hydro- 
dynamic interactions that shape the dispersion curve in 
the long- wavelength limit. In contrast, the short-ranged 
"collisions" between the droplets result in the opposite 
effect: at high q the density- wave packets propagate up- 
stream (u gx < 0, for qRd > 1). 

CONCLUSION 

In conclusion, we shed light on the spatiotemporal fluc- 
tuations generically observed in flow-driven confined sus- 
pensions. We devised a model microfluidic experiment 
which made it possible to track the individual positions 
of hundred of thousands of identical droplets, interact- 
ing hydrodynamically, in a shallow channel. We demon- 
strated that density fluctuations freely propagate at all 
scales and along all directions in the channel, irrespec- 
tive of the area fraction occupied by the droplets. In- 
troducing a kinetic theory we elucidated the two phys- 
ical ingredients that shape the dispersion curve of the 
linear excitations. At long wavelengths, non-dispersive 
sound modes propagate as a result of the dipolar hydro- 
dynamic interactions caused by the motion of the con- 
fined droplets. Conversely, for wavelengths smaller than 
ten particle diameters the strongly dispersive nature of 
the density waves stems from the combination of both the 
excluded volume interactions and the hydrodynamic cou- 
pling. A challenging perspective to this study would be 
to consider the case of self-driven particle such as rigidly 
confined motile cells, and chemically propelled colloids 
and droplets which also display intriguing coherent struc- 
tures at all scales [9~1|2HI22]. 

This work was partly funded by the Paris Emergence 
research program, and C'Nano He de France (DB). We 
acknowledge stimulating interactions with C. Savoie and 
M. Guerard. We thank Bertrand Levache for help with 
the experiments. 

MATERIALS AND METHODS 

Microfluidics and image acquisition 
The microfluidic device is a double etched fused sil- 
ica/quartz custom-design chip (Micronit Microfluidics). 
It consists of a conventional flow-focusing junction fol- 
lowed by a dilution module and a large channel. We 
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used deionized distilled water with 0.1% SDS surfactant 
and fluorescein for the dispersed phase, and hexadecane 
(Sigma- Aldrich) for the continuous phase. The channel 
surface were cleaned in a UV/Ozone cleaner prior to the 
experiments, thereby making the glass surface highly hy- 
drophilic. The injection into the channels is controlled 
by high-precision syringe pumps (Nemesys low pressure 
dosing modules, Cetoni). The device was mounted on 
a Nikon AZ100 upright macroscope. A Basler Aviator 
av2300-25gm (4 Megapixel, 8bit) camera was used to 
record the movement of the droplets in a field of view 
of 2.70mm x 1.03mm at the center of the channel. The 
frame rate was set at 45Hz. 

Data analysis 

From our movies, we detected the positions of the 
droplets to sub-pixel accuracy using a built in routine in 
the image processing program Image J. The particle ve- 
locities were then computed using the Matlab version of 
the tracking software developed by David Grier, John C. 
Crocker and Eric R. Weeks [5D] . We used a custom Mat- 
lab routine to compute the power spectrum of the density 
fluctuations. First, we computed analytically the spatial 
Fourier transform of the density field at all time, using a 
Gaussian shape function of width i?d/15 (see main text), 
which typically corresponds to the precision on the parti- 
cle position detection. Then, the temporal Fourier trans- 
form was computed numerically. To obtain the dispersion 
relations curves, we detected for each pulsation, oj, the 
wave- vector corresponding to the maximum of the power 
spectrum at this given pulsation. 
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